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We assess the predictive accuracy of perturbation theory based estimates of changes in covalent 
bonding due to linear alchemical interpolations among molecules. We have investigated a bonding 
to hydrogen, as well as cr and tt bonding between main-group elements, occurring in small sets 
of iso-valence-electronic molecular species with elements drawn from second to fourth rows in the 
p-block of the periodic table. Numerical evidence suggests that hrst order estimates of covalent 
bonding potentials can achieve chemical accuracy if (i) the alchemical interpolation is vertical (fixed 
geometry), (ii) involves molecules containing elements in the third and fourth row of the periodic 
table, and (iii) a reference geometry is optimized. In this case, changes in the bonding potential 
become near-linear in coupling parameter, resulting in analytical predictions with very high accuracy 
(~1 kcal/mol). Second order estimates deteriorate the prediction. If initial and hnal molecules differ 
not only in composition but also in geometry, all estimates become substantially worse, with second 
order being slightly more accurate than first order. The independent particle approximation to 
the second order perturbation performs poorly when compared to the coupled perturbed or hnite 
difference approach. Taylor series expansions up to fourth order of the potential energy curve of 
highly symmetric systems indicate a finite radius of convergence, as illustrated for the alchemical 
stretching of H J. Numerical results are presented for (i) covalent bonds to hydrogen in 12 molecules 
with 8 valence electrons (CH 4 , NH 3 , HjO, HF, SiH 4 , PH 3 , HaS, HCl, GeH 4 , ASH 3 , HaSe, HBr); (ii) 
main-group single bonds in 9 molecules with 14 valence electrons (CH 3 F, CH 3 CI, CH 3 Br, SiHsF, 

SiHaCl, SiH 3 Br, GeH 3 F, GeH 3 Gl, GeHsBr); (iii) main-group double bonds in 9 molecules with 
12 valence electrons (CHaO, CHaS, CHaSe, SiHaO, SiHaS, SiHaSe, GeHaO, GeHaS, GeHaSe); (iv) 
main-group triple bonds in 9 molecules with 10 valence electrons (HCN, HOP, HCAs, HSiN, HSiP, 

HSiAs, HGeN, HGeP, HGeAs); (v) H/ single bond with 1 electron. 


I. INTRODUCTION 

Solving Schrodinger’s time independent equation for 
the unperturbed electronic ground-state within the Born- 
Oppenheimer approximation yields the potential energy 
surface (PES) of any molecule as a function of nu¬ 
clear charges {Zj} (stoichiometry), nuclear positions 
{R/} (geometry), and number of electrons N (molec¬ 
ular charge) .1^ The PES plays a fundamental role in 
chemistry and elsewhere because many properties can 
be derived from it. While one can study efficient ways 
of predicting the PES of single compound^^^ efficient 
estimates of PES of ensembles of molecules are more 
useful (and challenging) in the context of virtual com¬ 
pound design efforts.^h^ These efforts typically attempt 
to search chemi cal co mpound space (CCS) spanned by 
{{Zj}, {Rj}, for novel materials with desirable 

properties. As such, accurate yet efficient quantum me¬ 
chanics (QM) based PES estimates hold the key f or 
successful rational compound design applications 
While many inexpensive semi-empirical QM methods are 
available, for this study we restri ct ourselves to first prin¬ 
ciples in the spirit of Refs.^^^UlHlII jviQj-e specifically, we 
investigate the application of “alchemical” coupling to 
the problem of efficiently estimating the PES of new 
molecules using Taylor series expansions in CCS, rather 
than empiricism. 

The alchemical coupling approach can be related to 


grand-canonical ensemble theory (Widom insertion)P^lt^, 
and has been well established for empirical force-field 
based molecular dynamics studiesUsing QM, al¬ 
chemical changes are less common despite E. B. Wil¬ 
son’s early proposal of variable Z, back in 1962.1^ Within 
QM, any two iso-electronic molecules in CCS can be cou¬ 
pled “alchemically” through interpolation of their exter¬ 
nal potentials. Here, we have investigated if alchemical 
predictions can be used to model the PES of covalent 
bonds occurring in small closed-shell molecules made up 
from main group elements. We have limited ourselves 
to covalent hydrogen bonds, as well as single, double, 
and triple bonds in molecules with no more than 14 va¬ 
lence electrons. We present and discuss numerical ev¬ 
idence for the following set of observations: First or¬ 
der Taylor-expansions of covalent bonding potentials can 
reach chemical accuracy (~1 kcal/mol) if two conditions 
are met. Firstly, the alchemical change has to be “ver¬ 
tical”, meaning that initial reference molecule as well as 
final target molecule have to possess the same number of 
atoms located at the exact same positions. Secondly, all 
elements involved in the alchemical change, i.e. all {Zj} 
destined to vary, have to occur late in the periodic table. 
Second order Taylor-expansion based predictions are less 
accurate than first order predictions if these conditions 
are met. If reference and target molecule have different 
geometries, the predictive power of the first order Taylor 
expansion substantially deteriorates, while second order 
estimates based on coupled perturbed Kohn-Sham equa- 
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tions offer some improvement, however, without reaching 
chemical accuracy. Second order estimates based on the 
independent particle approximation result in Taylor ex¬ 
pansion estimates that are even worse than first order 
estimates. For highly symmetrical alchemical changes, 
such as the dissociation of HJ, a finite radius of conver¬ 
gence is found. 

In Sec. [IT] we briefly summarize the framework of al¬ 
chemical derivatives within Hartree-Fock and density 
functional theory (DFT) as well as our notations. Nu¬ 
merical estimations of covalent bond stretching of small 
molecules are presented and discussed in Sec. |IIII Ex- 
tending previous work on alchemical perturbation^^oOll 
we discuss alchemical energy derivatives with respect to 
vertical transmutation, interpolating only the identity of 
the atoms while keeping the geometry fixed. Estimates 
of single, double and triple bonds are included as an ap¬ 
plication. We also report numerical results for alchemical 
stretching of chemical bonds using non-vertical transmu¬ 


tations. Finally, conclusions are drawn in Sec. IV 


be scaled down to/up from zero if the number of atoms 
in one molecule is smaller. Under isoelectronic condi¬ 
tions, the A-dependent terms in the coupling Hamiltonian 
(Eq. ([^) are the electron-nucleus and nucleus-nucleus in¬ 
teraction operators. 
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Since different pairing schemes result in different v\{r) 
and Va, it is obvious that the alchemical perturbation 
is alignment dependent. To investigate the behaviour of 
higher order corrections and the effects of varying geom¬ 
etry/stoichiometry, we neglect all relaxation effects for 
vertical iso-valence-electronic changes (see Sec. IIG). 


B. First order derivative 


II. METHOD 


A. Taylor expansion in CCS 


A Taylor expansion in CCS can be constructed 
with the exclusive knowledge acquired by solving 
Schrodinger’s equation for some reference molecule, with 
Hamiltonian H^, 


E{AX) = E^ + AXdxEx 


A=0 


AA2 

2 


diEx 



Derivatives of the total potential energy can be ob¬ 
tained by coupling a reference Hamiltonian to some tar¬ 
get Hamiltonian, Hj', such that Hx transforms into 
Hj' 


Hx = {l- X)H^ + XHt, (2) 

as the coupling parameter A goes from 0 to 1. And conse¬ 
quently, d^Ex = d'^{Hx), with dxHx = Ht - = H' 

being the alchemical perturbing Hamiltonian. If these 
derivatives can be computed, E^ can be estimated ac¬ 
cording to Eq. Q by setting AA = 1. Note that we cou¬ 
ple reference and target systems in a linear and global 
fashion. This is an arbitrary choice, non-linear and lo¬ 
cal interpolation functions could have been chosen just 
as well. In fact, in Ref. Eol an empirical quadratic in¬ 
terpolation function is found to yield superior results for 
first order predictions of highest occupied molecular or¬ 
bital (HOMO) eigenvalues. In this study of alchemical 
changes of covalent bonding, we begin with linear and 
global interpolations, future work might deal with alter¬ 
native functions. 

Given a pair of isoelectronic reference/target systems, 
described by {{Zf}, {R^'}, N} and {{Zj}, {Rj}, N} re¬ 
spectively, one can couple the two systems such that cer¬ 
tain Zf^ and Zj are paired. Note that Zf or ZJ can 


The first order derivative of the energy with respect 
to an alchemical interpolation parameter connecting any 
two iso-electronic molecules, can be computed accord¬ 
ing to the Hellmann-Feynman theorem,!^ as shown for 
molecular HOMO eigenvalues ,1^ 

dxEx = {dxHx)x = j px{r)dxvxir) + dxVx, (4) 

where PA(r) denotes the electron density, dependent on A. 
At A = 0 we have PA(r) = PR(r), which is independent 
of the target system. As such, the first order deriva¬ 
tive can be calculated with a single reference density 
and without additional self-consistent field (SCF) calcu¬ 
lation for any target system. In several circumstances, 
Taylor expansion estimates using first order alchemical 
derivatives have shown good accurac y for the rap id pre¬ 
diction of properties throughout CCSP^IlilSlISlMl ggjj_ 
eral, however, first order derivatives might not be suffi¬ 
cient. Taking higher order derivatives into account might 
offer higher accuracy, assuming Eq. 0 converges rapidly. 


C. Second order derivative 

Differentiation of Eq. 0. based on linear interpolated 
Hamiltonian in Eq. ([^, yields 

dlEx = J dr (dxp(r))(dxvx(r)), (5) 

requiring the density response due to the alchemical per¬ 
turbation. Again, at A = 0 this amounts to the density 
response of the reference system. Evaluation of Eq. 0 
implies a differing density response for each target sys¬ 
tem. We have considered three approximations to dxp 
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including second order perturbation theory with inde¬ 
pendent particle appro ximati orP^ (IPA), coupled per¬ 
turbed (CP) approachesas well as finite difference 
approximation (FD). Note that Eq. (§ can be rewrit¬ 
ten as dlE = J drdr'{dxv{r)){dxv{r'))-g:^^^££^, where 

Sv{r)^(T') ~ is the static linear response func¬ 

tion or susceptibility, well established within conceptual 
QpqUlHll] 

Perturbation theo ry pro vides ways to estimate 
9APA(r)PSl Within iPAPsESIlZl^ the static density response 
for a close-shell system is approximated by 


D. Higher order derivatives 

Mpller-Plesset (MP) perturbation theorj^SlMI is used 
to estimate correlation energy corrections based on con¬ 
verged Hartree-Fock results. The derivation of higher 
order corrections in MP theory are equivalent to the 
order alchemical derivative. Here, instead of the two- 
particle operator for electron-electron interaction as per¬ 
turbation in MP theory, the alchemical perturbation op¬ 
erator Ht — Hfi can be used. Within IPA, the MP for¬ 
mula can be directly applied to obtain any order 
derivative. 


dxpx{r) 
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E. Predicting changes in covalent bonds 


where denote the occupied molecular or¬ 

bitals (MOs) and their eigenvalues, while {^ai^a} de¬ 
note the unoccupied counterparts. IPA neglects the 
influence of the alchemical perturbation on th e Hartree 
and exchange-correlation (xc) potentials.^^^^^ Note that 
Eq. Q becomes numerically exact for 1-electron system 
with converged basis set within Hartree-Fock approxima¬ 
tion, because of the absence of Coulomb and xc interac¬ 
tion between electrons. 

Recently, Yang, Cohen, De Proft and Geerlings derived 
an expression of the density response that also includes 
the dependence of Coulomb and xc potential,!^ the CP 
approach,!^ 

dxpx{r) = -4 

ij ah 

x{M~^)ia,jb J dr' (j)jiv')(j)b(r')dxvx{r'), 

( 7 ) 

where the matrix elements of M are 


For the study of covalent bonds we focus on the changes 
in binding potential due to alchemical coupling. We con¬ 
sider the difference in total potential energy between two 
bounded atoms at two arbitrary interatomic distances d 
and do, 

AE{d,do) = E{d) - E{do). (9) 

If, for example, do is large and d is the geometry mini¬ 
mum, AE becomes the bond dissociation energy. We are 
interested in changes of AE{d, do) as a function of d due 
to alchemical changes for a fixed do. More specifically, we 
couple a reference to target system via the corresponding 
Hamiltonians yielding expectation values as a function of 

A, 

AEx{d,do) = Exid)-Ex{do) (10) 

= (Huid) + XiHrid) - Huid))) 

— (F7fl(do) + A(idT(do) — Hii{do))). 
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In the limit of Jiajb 0 and Xia jb 0, Eq. (|^ and 
Eq. 0 are equivalent. 

Alternatively, one can also introduce an explicit small 
perturbation and converge the new density at A A ^ 1. 
The density response can then be estimated via ED, 
dxpir) « practice, instead of starting the 

SCF for the perturbed system from atom based initial 
guesses, we restart with PR(r) resulting in convergence 
within few SCF steps. 


As A goes from zero to one, the two components in 
Eq. 0 change from reference {EB,{d), E^{do)) to tar¬ 
get {E'T{d),ET{do)) compound. The truncated Taylor 
expansion based estimate of the target compound’s po¬ 
tential is then obtained via, 

AET{d,do) - A4'")(d,do) = AER{d,do) (11) 

m - 

+ EM^AAAA(d,do), 

k=l 


where the superscript m stands for Taylor expansion with 
m terms, as a function of bond-length d for vertical al¬ 
chemical changes. Since AEt is the property of interest, 
the subscript T, A, and the dependency of do will be 
omitted for the rest of this work, unless otherwise noted. 
In this study we investigated orders up to m = 4 for the 
stretching of HJ, and up to m = 2 for all other molecules. 
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For a fixed do, first and second order estimation are 

= (Enid) + d^E^id)) (12) 

— {Enido) + d\E\{do)), 

AF;(2)(d) = (Enid) + d^Exid) + ^dlE^id)) (13) 

— {En{do) + d\E\(do) + -dfE\{do)). 

Since d and do in Eq. are arbitrary, one can infer the 
binding curve via scanning d for any fixed do- The pre¬ 
dictive power, however, happens to dependent on do- For 
this reason, we optimize do such that the integrated error 
in dissociation region is minimal. As shown in Fig. an 
empirical linear relationship exists between equilibrium 
bond length of target molecule and dopt 


dopt«0.76dJq + 0.97A. (14) 


do is determined according to Eq. (14) for all verti¬ 
cal changes. If djq is not known it can easily be esti¬ 
mated with semi-empirical quantum chemistry methods. 
For non-vertical changes, we fix do = deq to the equi¬ 
librium distance in th e refere nce molecule, resulting in 
AF;('")(deq) = 0. Eqs. become 


AE^^\d) = (Enid) + dxE^id)) (15) 

-(ER(deq) -f d),Ex{deq)), 

A£;(2)(d) = (Enid) + dxE^id) + ^dlEx{d)) 

~(d?R(deq) + d\E\{deq) + -O^Ex^deq)) ■ 


F. Error measures 


been set to correspond roughly to the inflection point, 
due to the issues of a single reference such as DFT 
method for describing covalent bond-dissociation. This 
shortcoming is also evident from comparison of DFT to 
CCSD(T) curves shown in Fig. 

Note that this aspect is irrelevant for the alchemical 
predictions: If a more reliable reference method had been 
used the error integration could easily be expanded to 
include the entire dissociative tail. These four quanti¬ 
ties provide a numerical indication of how good a pre¬ 
diction is. For a perfect prediction one would expect 
(AEeq, Adeq, Aw,IE) = (0,0,0,0). Note that we com¬ 
pare the predictions to DFT. This is an arbitrary choice, 
any other QM method could have been applied just as 
well. 


G. Computational details 

Alchemical interpolations of molecules containing ele¬ 
ments from different rows in the periodic table can still be 
iso-electronic if effective core or pseudo potentials (PPs) 
are used, resulting in a constant number of valence elec- 
tronJi^. For example, one can couple carbon to silicon us¬ 
ing just four valence electrons. Non-local PPs are widely 
used to mimic the presence of core electrons in atomPSl^ 
and are amenable to the tuning of a wide range of proper¬ 
ties includiM dispersion forces, band-gap, or vibrational 
frequencied^i**^. The non-local external potential nA(r) 
in Eq. ([^ then becomes 

Ni 

UA(r,r') = X! “ A)?;f(r,r') -f AuJ(r,r')), (17) 

I 


Eor bond lengths, we quantify the predictive power of 
the Taylor expansions by evaluating the deviation of pre¬ 
diction from the DET bond length Adeq = — deq, 

where di^^ stands for the predicted equilibrium distance 
of AE^"^\ We calculate the deviation of the predicted 
energy at di™^ from the DET energy at the DET mini¬ 
mum, AEeq = AA(™^(di™^) — AE{deq). The deviation 
in harmonic vibration frequency. Aw = w*-™^ — w, of 
the stretching bond is also included in order to quan¬ 
tify the accuracy of the stiffness of the predicted binding 
potential. The vibration frequency is computed from the 
curvature of cubic spline interpolated binding potential, 

w = where k^q = 9^AA(deq) and fi is the re¬ 

duced mass. Finally, we measure the integrated error 
(IE) for the dissociative tail, defined as 


IE = 


Id 


max 


1 




dx\AE‘^"^\x) - AE{x)\, (16) 


for vertical iso-valence-electronic changes. Note that 
while in principle one would like dmax —t oo, dmax lias 


where vf and vj are PPs for Zf- and Zj respectively. 
Note that t’A(r, r') in Eq. 0 and vx (r) in Eq. ([^ result 
in different coupling Hamiltonians, and therefore differ¬ 
ent A-dependencies of the energy and its derivatives. 

All results have been obtained within the Born- 
Oppenheimer approximation, where nuclei are clamped, 
nuclear repulsion Va is decoupled from the electronic 
wavefunction, and is added as a geometry- and A- 
dependent constant to the electronic energy. Nuclear- 
nuclear repulsion energy is computed automatically by 
most QM codes. However, it must be removed and 
recomputed independently for Va according to Eq. 
to avoid self-repulsion between transmutating atoms. 
Throughout the present study, standard atomic and 
plane-wave basis functions, linearly interpolated PPs, as 
well as the PBE xc potential^ within KS-DET is used. 
The scanning of O.SA < d < S.OA is carried out with 
increments Ad = 0.1 A. Eor each prediction order m, 
AE^™^ (d) are interpolated with cubic splines, from which 

the stiffness 9^Ai?^"*^(de™^) = k^q is computed. All den¬ 
sity volumetric data is printed into Gaussian CUBE files, 
from which integrated density slices are calculated. 
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1. Details for vertical iso-valence-electronic changes 


Numerical results for vertical iso-valence-electronic al¬ 
chemical changes (discussed in Sec. Ill A and IIIBl have 
been obtained with cpmeP, a plane wave basis with 100 
Ry cutoff, and Goedecker PPs.l^^^The periodic super¬ 
cell size is 20 x15x 15A^, and one heavy atom is fixed 
at (7.5 A, 7.5 A, 7.5 A) while the stretching atom shifts 
along -l-x-axis. For each geometry, heavy atoms are mu¬ 
tated to other elements in the same column of the peri¬ 
odic table while all H are fixed at the same location as 
in the reference compound. 


Since Eq. 0 is a non-local operator, Eqs. Q and 
([^ need to be converted to wavefunction expressions. 
The hrst order derivative for the Hamiltonian i7u_>T is 
evaluated using RESTART files in which the reference 
compound’s density and wavefunctions have been stored: 
dxE = (9a77)r = Et[pyi] - Er[pr]. And the second 
order derivative is evaluated correspondingly relying on 
ED, dlE « ^ ^ Q Q 5 ^ Wavefunc¬ 
tions of reference compound are used for while 

AEpp is evaluated by FD with linearly interpolated PPs 
parameter. 


Coupled-cluster results (CCSD(T)) obtained for HCl 
and HBr have been computed using GaussianO^^ in 
aug-cc-pVTZ-l^ basis, and default input parameters. 


2. Details for non-vertical iso-electronic changes 


Numerical results for non-vertical iso-electronic al¬ 
chemical changes have been obtained using atom cen¬ 
tered basis-sets. Restricted open-shell Hartree-Fock cal¬ 
culations have been carried out using Cartesi an aug - 
cc-pVTZ basis seiP^ for (Discussed in Sec. IIIDll. 
Eq. (§ and higher order derivatives are evaluated ana¬ 
lytically by Gaussian expansion of MOs. Reference ge¬ 
ometry is first relaxed by GaussiamO^SH ^nd the con¬ 
verged MO coefhcients are extracted to evaluate orbital 
integrals. NWCheiiP3 is used to scan AE as a function of 
A in Fig. [^d) along alchemical path with discretization 
AA = 0.01. It is done by reassigning nuclear charges in 
the system. 

Non-vertical alchemical changes in 10-electron 
molecules (discussed in Sec. HID 2) have been calculated 
using the uncontracted Cartesian Def2TZVP basis 
setP^. Uncontracted neon basis is used for second row 
heavy atoms. Additional hydrogen basis functions are 
placed along the stretching pathway, from d = 0.5 A 
to d = 3.0 A in increments Ad = 0.1 A. All systems 
with integer nuclear charges have been calculated using 
GaussianOd^ while systems with fractional nuclear 
charges have been calculated using NWCheiJ^ with 
discretization A A = 0.01. For each 0 < A < 1, the 
atomic density for SCF initial guess iterates through {C, 
N, O, F, Ne} to ensure convergence. In all Gaussicin and 
NWChem calculations we used Cartesian/Real spherical 


harmonic basis functions. 



d [A] 


FIG. 1. AE is shown as a function of d in Eq. (111. White 
background panels: True (black circles), first (red squares) 
and second (blue triangles) order predictions of changes in 
the covalent bond of hydrogen due to vertical alchemical in¬ 
terpolations. Gray background panels: The true potentials of 
the reference compounds employed for the predictions for the 
first and second order predictions. 
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III. RESULTS AND DISCUSSIONS 


A. Vertical iso-valence-electronic changes of X-H 

1. Predicted potentials 

Using Taylor expansions binding potentials have been 
estimated for covalent bonds involving hydrogen (X-H) 
for the following 12 molecules with 8 valence electrons: 
CH 4 , NH 3 , H 2 O, HF (second period); SiH 4 , PH 3 , H 2 S, 
HCl (third period); and GeH 4 , AsHs, H 2 Se, HBr (fourth 
period). Numerical results for vertical first (red) and sec¬ 
ond (blue) order truncated Taylor series estimates feature 
in Fig. They measure the change in X-H binding en¬ 
ergy as one goes from reference to target compound. 

We first note that the entire potential is repro¬ 
duced in semi-quantitative fashion for all combinations 
of reference/target molecules. The precise predictive 
power strongly depends on the choice of reference/target 
molecule pair, on the choice of doj and on the expansion 
going up to first or second order. First order estimates 
among molecules with elements from third or fourth row 
are very accurate (See Fig.[^ bottom and mid row in mid 
and bottom panel, respectively). By contrast, predict¬ 
ing, or starting with, second row elements consistently 
yields worse results. Inclusion of second order correc¬ 
tions does not necessarily lead to improved performance. 
Second order truncated Taylor series estimates only yield 
more accurate predictions than first order when the refer¬ 
ence molecule contains heavier elements than the target 
molecule. For example, if we predict HF using HBr as 
a reference, the second order prediction is more accurate 
than first order. Making the inverse prediction (i.e. HBr 
from HF), however, first order is more accurate than sec¬ 
ond order. 

The performance of truncated Taylor series dramati¬ 
cally varies depending on the choice of the do value. The 
top panel in Fig. [^illustrates this for A£’(2A,do) for 
HF—;>HBr as a function of A, once with dg = 0.94A— 
the equilibrium bond length of HF— and once with 
do = I. 57 A, a value for do which happens to linearize AU 
in A. While the coupling path of total energies is hardly 
distinguishable for £'(2A), i?(1.57A), and i?(0.94A), AU 
is strongly dependent on the choice of dp. By choosing 
do = I.57A, Ai?(2A, I.57A) in Eq. (10 1 becomes nearly 
linear, while plotting AE(2A,0.94A) reveals substantial 
curbing. This is why, when choosing the right do, first 
order predictions of AU can be very predictive. 

The top panel in Fig. also explains why second or¬ 
der estimates can be worse than first order, and why 
this changes for the reverse coupling: On the side of the 
lighter element (A —> 0), a weak convexity is noticable 
in AE (blue line), despite the overall concavity of the 
path. The presence of inflection points will always lead 
to a deterioration of second order predictions, resulting 
in a more accurate first order estimate. On the other 
side (A —)■ 1), no such inflection point exists and the sec¬ 
ond order term results in the expected improvement of 



FIG. 2. Alchemical coupling of HF (A = 0) to HBr 
(A = 1). TOP panel: E and AE where deq = 0.94 A (red) 
denotes the equilibrium bond length of reference molecule 
HF, and dopt = 1.57 A (blue) linearizes AE. BOTTOM 
panel: Integrated valence electron density difference slices 
between H-X at d = 2A and at dopt and deq, respectively, 
APx{x) = / dydz[p\{v,d) — px{r,do)]. Dependence on A is 
shown for the same vertical interpolations, Left and Right cor¬ 
responding to the non-linear (red) and linearized (blue) AE 
curves in right-hand TOP panel. Density changes at heavy 
atom and hydrogen positions are highlighted as red dashed 
and dotted/dash-dotted lines, respectively. 


the prediction. For the other compound pairs shown in 
Fig.[T] similar observations can be made for the attrac¬ 
tive part of the bonding potential (see also AU(2A, do) 
curves in supplementary materials). We have found that 
the inflection point near A = 0 occurs only when the ref¬ 
erence molecule has the lighter element. Conversely, no 
inflection point has been observed for atoms transmutat- 
ing upward the column. We believe that this behaviour is 
due to the specifics of the employed PPs. Future studies 
will show why this is the case, and if similar trends hold 
for other PPs. 


2. Integrated error 


Prediction errors for energy minima, equilibrium bond 
lengths, force constants, and integrated error in d issoci- 
ation region (calculated as described in Sec. HF) have 
been obtained for all predictions in Fig. and are listed 
in Table. |T] The results lend quantitative support for the 
observations articulated above. In particular, the results 
suggest that chemical accuracy can be obtained when 
using first order Taylor series based estimates among 
compounds containing third and fourth row elements. 
Second order based predictions are always worse except 
when a molecule with heavier element is used as a refer¬ 
ence to predict a molecule with lighter one, for example 
HBr^HF. 
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The best prediction performance is found for first or¬ 
der based estimates using reference molecules containing 
third row elements (riR = 3) in order to predict target 
molecules made up of fourth row elements (riT = 4). 
The overall average deviation from reference bonding po¬ 
tential energies and integrated error are ~2.5 kcal/mol. 
Corresponding predictions of equilibrium distances devi¬ 
ate at most by O.OSA, and the vibration frequencies de¬ 
viate no more than 32 cm“^. Second order estimates for 
the same third and fourth row combinations give slightly 
worse results. The worst predictions are found if the 


coupled molecules skip a row, i.e. involve elements from 
second and fourth row—for first as well as second order 
truncated Taylor expansions. This is not surprising as 
the central atom’s electron density must accommodate 
the most severe contractions/expansions for such inter¬ 
polations. Moreover, second order overcorrections can 
also be found (Table. |l]) whenever molecules containing 
fourth row elements are used to predict molecules con¬ 
taining third row elements: Both, predicted energy min¬ 
imum and equilibrium bond length, show negative devi¬ 
ations. 


TABLE I. Error measures for first (left) and second (right) order predictions of vertical iso-valence-electronic alchemical 
changes of covalent bond potentials in X-H—^A-H. The compound pairs are arranged in the same way as in Fig. i.e. each 
box corresponds to an unshaded panel containing a predicted potential. Unit are [kcal/mol] for AEeq and IE, [A] for Adeq, 
and Alo in [cm“^]. Avg corresponds to averaged signed error for each row with corresponding unit. Periods of X and A are 
denoted by their respective primary quantum number nR,nT = {2,3,4}. IV, V, VI, VII represent the columns of X and A in 
the periodic table. 

nn = 2 nn = 2 


riT 

IV. 

V. 

VI. 

VII. 

Avg 

A£;eq 
Q Arfeq 

= 3 Ac. 

IE 

10.55 

0.12 

-46.9 

13.01 

17.19 

0.14 

-109.6 

20.00 

22.16 

0.15 

-188.5 

27.43 

24.05 

0.15 

-213.7 

30.18 

18.49 

0.14 

-139.7 

22.65 

AEeq 

A Adeq 

= ^ Ac 

IE 

10.76 

0.15 

-22.8 

13.19 

19.97 

0.20 

-137.8 

28.19 

23.77 

0.24 

-164.1 

36.84 

25.74 

0.21 

-196.8 

39.42 

20.06 

0.20 

-130.3 

29.41 



AE ^^^, tt-r = 4 


A.E'eq 
r> Adeq 

IE 

38.08 

0.26 

-331.4 

60.39 

44.22 

0.31 

-455.1 

78.48 

54.03 

0.33 

-636.6 

106.16 

61.07 

0.33 

-738.0 

121.95 

49.35 

0.31 

-540.3 

91.75 

A.E'eq 

1.26 

2.43 

3.81 

5.39 

3.22 


0.01 

0.02 

0.02 

0.04 

0.02 


-14.4 

-15.0 

-31.0 

-60.4 

-30.2 

IE 

1.16 

2.31 

3.56 

5.13 

3.04 


3. Alchemical predictions do not commute 

We note the asymmetry in the predictive power of first 
order based predictions which is due to the lack of com¬ 
mutation: In general i9a^'U=o 7 ^ d\E\x^i, except if ref¬ 
erence and target Hamiltonian happen to differ only by 


Tlx 

IV. 

V. 

VI. 

VII. 

Avg 

AEeq 

18.92 

26.34 

31.39 

33.94 

27.65 


0.25 

0.27 

0.28 

0.26 

0.26 

= 3 Ac" 

-53.6 

-144.1 

-233.7 

-256.4 

-172.0 

IE 

32.34 

41.13 

52.19 

56.39 

45.51 

AEeq 

18.14 

27.34 

30.60 

33.66 

27.43 


0.28 

0.34 

0.36 

0.34 

0.33 

= 4 Ac" 

-39.1 

-136.4 

-141.2 

-182.0 

-124.6 

IE 

32.30 

55.27 

67.82 

68.92 

56.08 


AE^^\ riR = 3 


AEe^ 
n Ac?eq 

^^ = 2 Ac 

IE 

12.33 

0.07 

-143.1 

12.95 

15.46 

0.08 

-202.4 

16.79 

19.08 

0.09 

-259.0 

21.26 

20.63 

0.08 

-312.3 

22.16 

16.88 

0.08 

-229.2 

18.29 

AEeq 

4.66 

7.24 

9.33 

10.90 

8.03 


0.04 

0.08 

0.09 

0.09 

0.08 

^^ = 4 Ac" 

-25.1 

-48.3 

-56.5 

-68.3 

-49.6 

IE 

4.83 

8.69 

10.69 

11.89 

9.03 


AE ^'^'^, riR = 4 


AEeq 
„ Ac/eq 
= 2 Ac 

IE 

19.91 

0.12 

-211.2 

22.83 

22.19 

0.12 

-276.1 

25.10 

30.44 

0.16 

-418.1 

38.88 

33.75 

0.15 

-487.4 

42.45 

26.57 

0.14 

-348.2 

32.31 

AEeq 

-1.30 

-3.64 

-5.32 

-6.49 

-4.19 


-0.01 

-0.03 

-0.03 

-0.04 

-0.03 

^^ = 3 Ac" 

12.3 

12.4 

24.0 

51.1 

24.9 

IE 

1.43 

3.07 

4.39 

5.23 

3.53 


translation, rotation, or parity (enantiomers, i.e. without 
accounting for parity violation). Within our restricted 
case of linear interpolations of iso-electronic systems, the 
perturbing potential does differ only by sign. The inte¬ 
gral over its product with the electron density, however, 
differs in general, i.e. (TJa — 77b )b = J dr pb{va — vb) 


































































f dr pa{va — vb) = —{Hb — Ha) A- As such, the error in 
estimating A based on B will not be the same as the error 
in estimating B based on A. Results in Table. |T] suggest 
that predictions downward the columns in the periodic 
table are more accurate than upward. For example, pre¬ 
dicting HBr using HF as a reference, a better estimate 
is obtained (error = -1-25.7 kcal/mol) than for predicting 
HF using HBr as a reference (error = -1-61.1 kcal/mol). 
Correspondingly, predicting HCl using HBr has an er¬ 
ror = -1-5.4 kcal/mol, while the prediction of HBr using 
HCl has only an error of -1-3.6 kcal/mol. Similar obser¬ 
vations hold for bond lengths, and force constants. The 
asymmetry is also illustrated in Fig. AE{d,do) is not 
necessarily symmetric with respect to A = 0.5 for a given 
choice of (d, do)- Consequently, truncated Taylor series 
based predictions from either end will not be equally ac¬ 
curate. 



d[A] 


FIG. 3. Alchemical predictions can be more accurate than 
approximated density functionals. Covalent binding poten¬ 
tials obtained from alchemical PBEO estimate (red squares) 
and ordinary PBE (green diamonds) for HBr (full symbols) 
and HCl (empty symbols). The alchemical estimate corre¬ 
sponds to Eq. ( [T^ using PBEO density of HBr (HCl) in or¬ 
der to predict HCl (HBr). For comparison, corresponding 
CCSD(T) results are shown as well (dashed). 


4- Chemical accuracy 


B. Vertical iso-valence-electronic changes involving 
single, double, and triple bonds 


We have seen that very accurate, yet inexpensive, first 
order alchemical estimates can be made for vertical al¬ 
chemical changes between third and fourth row elements 
according to Eq. Q— once the density is converged for 
a given reference molecule. Then, an interesting ques¬ 
tion is if the alchemical accuracy is on the same order of 
magnitude as common approximations made when solv¬ 
ing Schrodinger’s equation. We have investigated this 
point for alchemical coupling of HBr and HCl using hy¬ 
brid and generalized gradient approximated DFT. When 
using PBElP^ as the method for the reference compound, 
we find the first order based alchemical predictions ac¬ 


cording to Eq. (12) to be in better agreement with the 
PBEO results for the target compound than true gener¬ 
alized gradient based approximation PBEP^. Eig. il¬ 
lustrates this point for the covalent binding potentials of 
HCl and HBr calculated using PBEO, PBEO based ver¬ 
tical first order alchemical predictions, and PBE. For all 
interatomic distances in the dissociative tail, the alchem¬ 
ical prediction (squares) is closer to PBEO (circles) than 
PBE (diamonds). For the repulsive part of the poten¬ 
tial, the alchemical prediction is substantially better than 
PBE for HBr, and slightly worse than PBE for HCl. For 
comparison we also included CCSD(T) results. These 
results amount to numerical evidence that the predic¬ 
tive power of vertical alchemical predictions can exceed 
the accuracy of common DFT approximations for third 
or forth row elements—if a sufficiently accurate electron 
density is provided for the reference compound. 


1. Predicted potentials 


Having discussed covalent bonds involving hydrogen, 
we now turn to single (XH 3 -Y), double (XH 2 =Y), and 
triple (HX#Y) bonds among p-block elements. Since 
third row elements can either be alchemically compressed 
to the corresponding second row (n = 2 ) element in 
the same column, or expanded to the fourth row (n = 
4) element, we chose third row (n = 3) based refer¬ 
ence systems for single, double, and triple bonds, namely 
SiHaCl, SiH 2 S, and HSiP. The resulting eight alchemical 
paths are combinations of changing the Si atom (Si—>C, 
Si—^Ge) or its binding partner (Cl—J-F, Cl—^Br, S—>-0, 
S—^Se, P—tN, P—J-As). In Figs. first and second order 
alchemical predictions are shown for the bonding poten¬ 
tial using vertical transmutations from the three refer¬ 
ence molecules. 

More specihcally, single bonds predictions have been 
investigated for making predictions using SiHaCl as a ref¬ 
erence compound for the eight following molecules with 
14 valence electrons: CH 3 F, CH 3 CI, CH 3 Br (nx = 2); 
SiHsF, SiH 3 Br (nx = 3); and GeH 3 F, GeH 3 Cl, and 
GeH 3 Br (nx = 4). For double bonds, we have considered 
predictions for the following eight unsaturated molecules 
12 valence electrons and using SiH 2 S as a reference com¬ 
pound: CH 2 O, CH 2 S, CH 2 Se (nx = 2); SiH 20 , SiH 2 Se 
(nx = 3); and GeH 20 , GeH 2 S, and GeH 2 Se (nx = 4). 
And finally for triple bonds, we have studied the following 
eight molecules with 10 valence electrons and using HSiP 
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as a reference compound: HCN, HCP, HCAs (nx = 2); 
HSiN, HSiAs (nx = 3); and HGeN, HGeP, and HGeAs 
(nx = 4). 
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FIG. 4. Alchemical predictions of single (top), donble (mid¬ 
dle), and triple (bottom) bond potentials. Curves are shown 
for eight target systems (specified as insets), iso-electronic 
with reference molecule SiHaCl (upper panel), SiHaS (middle 
panel), and HSiP (bottom panel) True (black circles), first 
(red squares) and second (blue triangles) order vertical al¬ 
chemical predictions of heavy atom bond dissociation curves. 


Numerical results in Fig. indicate qualitatively cor¬ 


rect behavior for all predictions. Regarding quantitative 
performance, the accuracy of the alchemical prediction of 
AE{d,do) exhibits similar behavior as the one discussed 
above in the case of vertical changes in the hydrogen 
containing single bond: First order predictions (red) sys¬ 
tematically achieve strong predictive power whenever the 
change involves the coupling of the third row element to 
a fourth row element. Gorresponding second order pre¬ 
dictions (blue) deteriorate the accuracy due to inflection 
points near A = 0. If the coupling involves one lighter 
element from the second row, the prediction is no longer 
quantitative. However, in these cases, second order pre¬ 
dictions provide a slightly superior prediction. If both 
atoms are simultaneously transmutated to lighter atoms 
from the second row, e.g. SiH 3 Gl—>'GH 3 F, second or¬ 
der estimates over correct (change of sign) the first order 
prediction. In the case of one element transmutating up¬ 
ward the column, the other downward, the second order 
estimate is hardly distinguishable from the first order es¬ 
timate. We believe that the reason for this is that the 
coupling to the lighter element on the one site in the 
molecule yields the concave behavior leading to an im¬ 
provement in the prediction, while the coupling to the 
heavier element on the other site in the molecule yields 
the convex behavior with the inflection point, leading to 
a deteriotation of the prediction. Effectively, these two 
effects cancel each other and result in the same predictive 
accuracy as the one obtained for the first order estimate. 
This rationalization rests on the assumption that the dis¬ 
cussion of Fig.j^can be applied also to linear combination 
of effects at different transmutating sites. 


2. Integrated errors 


Above observations are consistent with the quantita¬ 
tive integ rated prediction error measures (definitions in 
Sec. IIF) summarized in Table. [TT} All first order based 
predictions of target molecules implying a transmutation 
downward the periodic table (columns 4/3, 3/4, 4/4) ex¬ 
hibit chemical accuracy with at most 1.83 kcal/mol devia¬ 
tion in minimal energy (GeH 2 S), at most 0.04 A deviation 
in bond length (GeH 2 Se), at most -12.82 [cm“^] devia¬ 
tion in wavenumber (SiH 2 Se), and at most 1.56 kcal/mol 
in integrated energy (GeHAs). The best performance is 
achieved in the case of changing SiH 3 Cl—)'GeH 3 Br with 
energy error AE = 0.6 kcal/mol and integrated IE = 0.9 
kcal/mol. Corresponding predictions of equilibrium dis¬ 
tance deviates 0.03 A with vibration frequency deviate 
-l.I cm“^. First order predictions do not yield quanti¬ 
tative predictive power for changes involving lighter el¬ 
ements (columns 2/3, 3/2, 2/4, 2/2, 4/2). The worst 
predictions are found for the simultaneous coupling to 
two lighter elements (colum 2/2) with 76.29 kcal/mol, 
0.35 A, -568.76 cm“^, and 41.45 kcal/mol deviation in 
minimum energy, bond length, harmonic frequency, and 
integrated energy (HGN). 

While for all first order predictions all mutual devi- 
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FIG. 5. Scatter plot of optimized reference bond length dopt 
versus equilibrium bond length of target molecule, denoted 
by djlq. Linear regression gives dopt = 0.76 deq + 0.97 A with 
MAE=0.11 A and RMSE=0.15 A. First order (red empty 
squares) and second order (blue empty triangles) d^pt of co¬ 
valent X-H bond stretching, as well as first order (red filled 
squares) and second order (blue filled triangles) dopt of X-Y, 
X=Y, X#Y stretching are shown, where -, =, # stand for single, 
double, and triple bonds. Some of the alchemical paths are 
highlighted by black arrows. All numbers are given in Tables 
I and II in the supplementary material. 


ations exhibit the same sign, second order corrections 
introduce the sign changes in minimum energy and bond 
length alluded to before, namely both third row elements 
couple to lighter elements from the second row (colum 
2/2). Second order predictions for this column are even 
worse than the corresponding hrst order predictions. Sec- 


D. Non-vertical iso-electronic changes 


We are not aware of any mathematical limitation on 
how to construct alchemical coupling paths under iso- 
electronic condition. In addition to the investigation of 
predicted PES of iso-electronic compounds with the same 
geometry, as discussed in Sec. |III A| and Sec. |IIIB[ we 
have also investigated if one can use only one reference 
calculation in order to estimate the entire PES through 
“non-vertical” interpolations. In other words, we have 
also assessed the applicability of the Taylor expansion 
of Eq. Q to non-vertical changes for varying geometry 
and/or atom types and numbers between reference and 
target molecule. 


ond order predictions only improve first order predictions 
in the case of columns 2/3 and 3/2, alas, not to a degree 
considered satisfying. 

In summary, for changes corresponding to columns 
4/3, 3/4, 4/4 first order based estimates yield chemi¬ 
cal accuracy. For changes corresponding to columns 2/2, 
first order based estimates are inaccurate but still better 
than second order estimates. For changes correspond¬ 
ing to columns 2/4, and 4/2, first order based estimates 
are similar to second order estimates, yet both are inac¬ 
curate. For changes corresponding to columns 2/3 and 
3/2, second order based estimates are inaccurate but still 
better than first order estimates. 


C. Empirical d, 


'Opt 


Above observations have been made for optimized do¬ 


lt should be noted that the choice of dg in Eq. (10) is 


crucial for linearizing the property of interest in alchem¬ 
ical coupling parameter A, and hence essential for the 
performance of the perturbation based predictions. We 
have found that the error minimizing dopt has an ap¬ 
proximately linear dependence on the target molecule’s 
equilibrium bond length deq, no matter if the reference is 
the hydrogen containing single bond, or a single, double, 
or triple bond involving p-block elements from second, 
third, or fourth row. Furthermore, the linear relation¬ 
ship is preserved, independent of the fact if predictions 
are made with first or second order estimates. This rela¬ 
tionship is shown in Fig. The parameters of a linear 
regression are specihed as well. The outlier in Fig. [^at 
target dgq « 1.35A and dopt ~ 1-OA is due to the second 
order prediction of SiHsCl—^CHaF, i.e. for the above dis¬ 
cussed worst case scenario (column 2/2) where a strong 
overcorrection has been found. 


1. Alchemical stretching ofYl'l 


We now turn to the case of alchemical stretching of 
in order to understand the effect of varying geometry 
on alchemical predictions. Since Hartree-Fock is numer¬ 
ically exact for one-electron systems, we have employed 
an atomic basis set in an ’’all-electron” (no PPs) calcu¬ 
lation within the following alignment scheme: One pro¬ 
ton is centered at Ri = (0,0,0), the other is aligned 
along the -fee-axis. The reference system corresponds to 
at its equilibrium bond length. Stretching is accom¬ 
plished not by pulling the atoms apart but rather by si¬ 
multaneous annihilation and creation of nuclear charges 
at = (deq, 0,0) and Rj = (d, 0,0), respectively. Once 
the SCF is done for deq, the entire binding potential can 
be estimated up to to = 4 order, using Eq. (Ill, by scan¬ 
ning through various d, i.e. setting dg = deq. 

Results are shown in Fig. |^a). Due to the variational 
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TABLE II. Summary of error measure in Fig. for first (upper table) and second order (lower table) based predictions. The 
reference molecules are in the left hand column as Hb. = {SiHsCl, SiH 2 S, HSiP}. The primary quantum number of the heavy 
atoms X = {C, Si, Ge} and Y = {N, P, As, O, S, Se, F, Cl, Br} in target molecules XH 3 -Y, XH 2 =Y, and HX^Y, respectively, are 
specified using the corresponding principal quantum numbers nx and uy in each column. The right hand column (nx = d/nv 
= 4), for example, corresponds to predictions of molecules GeHaBr, GeH 2 Se, and HGeAs, respectively. Error measures of each 
panel are collected in corresponding cell with unit [kcal/mol] for Ai?eq and IE, [A] for Adeq, and Ao; is in [cm~^]. 


Hn 

2/3 

3/2 

2/4 

2/2 

4/3 

3/4 

4/2 

4/4 


A£^eq 

9.04 

18.25 

10.37 

20.37 

0.95 

0.80 

15.83 

0.60 

0 

CO 

AcZeq 

0.22 

0.27 

0.31 

0.29 

0.03 

0.03 

0.33 

0.03 

X 

AiU 

-91.83 

-169.02 

-106.27 

-237.85 

-12.05 

-5.02 

-148.62 

-1.09 

zn 

IE 

10.14 

16.82 

12.98 

14.85 

1.18 

0.98 

16.74 

0.88 


A£^eq 

21.69 

30.78 

22.87 

52.34 

1.83 

1.28 

26.77 

1.42 

m 

CN 

Ac^eq 

0.21 

0.25 

0.29 

0.32 

0.02 

0.03 

0.28 

0.04 

m 

AiU 

-191.24 

-241.71 

-201.69 

-436.06 

-9.88 

-12.82 

-218.43 

-10.84 

IE 

15.16 

19.47 

14.94 

33.25 

0.92 

0.92 

17.18 

0.98 


A£^eq 

31.35 

26.34 

30.42 

76.29 

0.79 

1.32 

23.91 

1.28 

X 

AcZeq 

0.21 

0.23 

0.25 

0.35 

0.01 

0.02 

0.24 

0.03 


AiU 

-240.70 

-229.35 

-250.08 

-568.76 

-1.77 

-9.15 

-192.45 

-2.34 

IE 

22.28 

19.42 

23.73 

41.45 

0.93 

1.22 

18.35 

1.56 


A£;(2) 



A_£/eq 

7.79 

11.84 

13.29 

-103.39 

3.00 

2.54 

17.60 

2.68 

0 

« 

Ac^eq 

0.13 

0.13 

0.39 

-0.18 

0.09 

0.10 

0.34 

0.13 

X 

Auj 

-91.30 

-111.67 

-147.58 

1912.55 

-22.96 

-18.27 

-176.49 

-13.20 


IE 

6.48 

8.47 

16.19 

32.26 

3.36 

3.61 

17.63 

3.83 


A_£/eq 

12.46 

17.88 

28.09 

-101.05 

5.18 

4.59 

29.36 

5.36 

m 

CN 

Ac^eq 

0.09 

0.11 

0.33 

0.06 

0.08 

0.09 

0.26 

0.13 


Auj 

-113.32 

-128.03 

-258.05 

1290.29 

-36.73 

-30.13 

-256.16 

-27.49 


IE 

6.84 

9.58 

15.34 

33.66 

2.47 

2.75 

14.82 

2.42 


A_£/eq 

12.04 

14.92 

30.85 

-50.92 

3.05 

4.22 

21.95 

4.93 

X 

Ac^eq 

0.07 

0.11 

0.22 

0.14 

0.04 

0.07 

0.18 

0.10 


Auj 

-94.83 

-118.62 

-256.64 

747.69 

-11.84 

-29.35 

-182.63 

-22.36 

IE 

8.56 

9.45 

21.74 

28.53 

2.96 

4.58 

14.80 

5.71 


principle for linearly coupled alchemical Hamiltonians ,1^ 
> AE for all interatomic distances. Inclusion of 
second order term improves upon the first order predic¬ 
tion, yielding a reasonable binding potential. However, 
when including third and fourth order the performance 
deteriorates again with oscillating behaviour for varying 
order (Fig. [^a) and inset of (b)), as AE^^^ overshoot 
and AE^^'> over corrects. Overall AE^‘^'> gives the best 
prediction. 

To explain the oscillating behaviour in Taylor expan¬ 
sion order, we investigate in more detail how the system 
responds to alchemical perturbation. When A increases 
gradually from 0 to 1, the nuclear charge decreases from 
1 to 0 at R^, while increasing from 0 to 1 at Rj. Using 
the alchemical derivatives at A = 0, truncated Taylor se¬ 
ries based estimates are plotted along the true energy in 
Fig. I^b) as a function of A at d = 3 A. AE^^\ AE^‘^\ 
AUW, and AF;(4) are linear, quadratic, third order, and 
fourth order polynomials, respectively. Clearly, the trun¬ 
cated Taylor series will fail to converge to AE at A = 1 
due to a sharp change of AE at A « 0.9. This implies 
a strong nonlinear electronic response occurring late in 
the alchemical coupling regime, resulting in the oscillat¬ 


ing behaviour of the predicted PES in Figs.j^a) and (b). 
Note that while the sign of error alternates, the magni¬ 
tude of error also increases as one increases the order. 
Similar behaviour can be observed for other values of d. 

The energy gain starting at A « 0.9 is due to a rapid 
rearrangement of electron density for A > 0.9. This is il¬ 
lustrated in Fig. I^c) where the integrated electron den¬ 
sity Px{x) is plotted as a function of both A and x at 
d = 3 A. Cohen and Mori-Sanchez already pointed out 
for H^ the dramatic changes in electronic structure for 
infinitesimallysmall changes in nuclear charges at infi¬ 
nite distance.l^ One would expect this effect to intensify 
as more basis functions are taken into account. This be¬ 
haviour can be seen in Fig. [^c). The locations of the 
proton at origin Ri, as well as the location of the an¬ 
nihilated proton at R^, and created proton at Rj, are 
indicated by red lines. 

Further analysis shows that for A > 0.5, both ground 
and first excited state orbitals are localized: The elec¬ 
tronic ground state is localized at Ri while the first ex¬ 
cited state is localized at Rj. At A « 0.9, the two eigen¬ 
values become degenerate, resulting in a rapid change of 
the ground state density in order to meet the non polar 
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FIG. 6 . order truncated Taylor series of are de¬ 
noted by in (a) as a function of d at A = 1 and in 

(b) as a function of A at d = 3 A. Inset shows the error of 
at A = 1 (c) Integrated density Px (x) = f dydz p\ (r) 
where r = (x,y,z), is presented as a function of both A and 
a: at d = 3 A where the integrated density values at nuclei 
locations highlighted by red lines at x = 0 A, a; = 1.1 A, and 
X = 3 A, while contour lines are draw at the bottom, (d) 
HOMO/LUMO levels, denoted by eh and el respectively, are 
plotted as a function of A at d = 3 A. 


FIG. 7. Energy difference AE, first order truncated Tay¬ 
lor series AE^^\ second order truncated Taylor series cal¬ 
culated by coupled perturbed AE^p, and second order trun¬ 
cated Taylor series calculated by independent particle approx¬ 
imation are plotted as black circles, red squares, blue 

open triangles, and blue filled triangles respectively. Gou- 
pling Hamiltonians are arranged as follow: (a) GH 4 —>■ GH 4 , 
(b) NH 3 ^CH 4 , (c) GH 4 ^NHg, and (d) NH 3 ^NH 3 . In¬ 
sets of (b) and (c) show the zoom-out energy scale for overall 
landscape. 


symmetry requirement of HJ, by taking a linear combi¬ 
nation of both ground and first excited state. Note that 
there is no orbital node at midpoint, indicating a true 
ground state for a dissociated molecule. The degener¬ 
acy occurs for the system with fractional nuclear charges 
at A « 0.9. The dramatic change in density stabilizes the 
system in A, giving rise to the sharp decrease in energy 
in Fig. §b) , as A increases from 0.8 to 1. Perturbation 
theory for degenerate cases might be necessary to prop¬ 
erly account for this case. The degeneracy of the ground 
state and first excited state is shown for the eigenvalue 
crossing in Fig. [^d): The eigenvalue of the (highest) oc¬ 
cupied molecular orbital (HOMO) and lowest unoccupied 
molecular orbital (LUMO) are plotted as a function of A. 
The degeneracy breaks when ground state and first ex¬ 
cited state switch order, which results in a delocalized 
ground state. By contrast, note that the eigenvalues will 
not cross each other if the stretching is carried out by 
moving in real space. 

Crossing of eigenvalue surfaces limits the radius of 
convergence of alchemical Taylor expansion series within 
electronic ground-state theories. As a result, the Taylor 
expansion for this system is not convergent at A = 1, sim¬ 
ilar to well known cases in Mpller-Ploesset theory.tS^H^ 
For asymmetric alchemical interpolations, as exemplified 
for the follow ing examp les in this study, as well as in pre¬ 
vious studies,^ ^ * ^^ * ^° * ^^ l the energy is typically smooth in 
all A values, and derivative based expansions are expected 
to converge. 


2. Non-vertical iso-electronic changes in ten electron 
systems 

In the final section of this paper, we consider alchemi¬ 
cal non-vertical changes of molecules with ten electrons. 
More specifically, we present numerical results of non¬ 
vertical iso-electronic changes involving bond stretching 
in second row systems {CH 4 , NH 3 , H 2 O, HF}, using all 
electron DFT. The example has indicated that non¬ 
vertical changes can profit from second order estimates. 
Since exact analytical expressions are not available for 
systems with so many electrons, and since no PPs are 
involved, we have relied on approximative second order 
expressions IPA and CP, rather than on finite difference 
expressions (see Methods section above). 

3. Predicted potentials 

Fig. [7] illustrates the prediction of R-H covalent bond 
potentials for CH 4 and NH 3 , predicted from alchemical 
derivatives using the electronic structure obtained by a 
single SCF. As a reference system we used once the re¬ 
laxed CH 4 system (panels (a), (c)), and once the relaxed 
NH 3 (panels (b), (d)) geometry. For the chemical com¬ 
position of Rr being the same as and only the bond 
being stretched (Fig. [^a), (d)), the first order estimate 
constitutes an upper bound, i. e. it always overshoots due 
to the concave behaviour of AE as a function of A, also 
on display in Fig. ib). When also changing the chem¬ 
ical compositions from CH 4 and NH 3 or vice versa, the 
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TABLE III. Prediction errors using first and second order 
based alchemical estimates for non-vertical changes in ten 
electron systems. Deviation of predicted energy minimum 
from actual AEeq [kcal/mol], corresponding bond length de¬ 
viation Adeq [A], and vibration frequency Acaeq. Numerical 
results from first order second order with independent 

particle approximation Alfjp^, and second order with cou¬ 
pled perturbed AEq^ truncated Taylor series are presented. 
Average is calculated for each row in Avg column. 


Hr 

CH 4 

NH 3 

H2O 

HF 

Avg 

AFeq 

- 0.22 

-2.30 

-9.63 

-21.99 

-8.54 

CH 4 Adeq 

- 0.01 

0.12 

0.22 

0.30 

0.15 

Atj 

6061. 

3545. 

2641. 

1874. 

3530. 

A£/eq 

-4.33 

-0.04 

-4.34 

-14.01 

-5.68 

NH 3 Adeq 

-0.13 

0.002 

0.10 

0.19 

0.04 

Ao; 

4218. 

4268. 

4622. 

2449. 

3889. 

A£/eq 

-32.78 

-5.97 

- 0.0000 

-4.88 

-10.91 

H2O Adeq 

-0.27 

- 0.11 

0.0005 

0.09 

-0.07 

Ao; 

4360. 

3933. 

3886. 

2568. 

3687. 

A-^eq 

- 120.8 

-38.63 

-7.95 

0.0003 

-41.85 

HF Adeq 

-0.44 

-0.25 

- 0.10 

-0.0008 

- 0.20 

Atj 

5772. 

4047. 

3982. 

3560. 

4340. 


AFeq 

-8.99 

-437.5 

-935.2 

-1356. 

-684.4 

CH 4 Ac/eq 

-0.31 

-0.41 

-0.40 

-0.39 

-0.38 

Atj 

2207. 

9625. 

12610. 

14600. 

9761. 

AFeq 

-113.4 

-4.10 

-673.2 

-1394. 

-546.1 

NH 3 Adeq 

0.29 

-0.27 

-0.45 

-0.43 

- 0.21 

Acj 

470.9 

1471. 

11200 . 

12660. 

6452. 

A£/eq 

-198.7 

-90.18 

-0.007 

-809.4 

-274.6 

H2O Adeq 

0.16 

0.21 

0.007 

-0.43 

- 0.01 

Acj 

2949. 

1124. 

-1484. 

9157. 

2937. 

A_E/eq 

- 212.6 

-140.4 

-64.95 

0.002 

-104.5 

HF Adeq 

0.05 

0.11 

0.15 

- 0.001 

0.08 

Acj 

6035. 

4110. 

1975. 

-1177. 

2736. 




AFeq 

-0.39 

-2.74 

-19.70 

-41.12 

-15.99 

CH4 Adeq 

- 0.02 

-0.07 

- 0.10 

- 0.10 

-0.07 

Atj 

1663. 

-495.7 

2116. 

3584. 

1717. 

AFeq 

0.90 

-0.04 

-1.84 

-15.17 

-4.04 

NH3 Adeq 

-0.04 

0.007 

-0.06 

- 0.10 

-0.05 

Alj 

2133. 

140.4 

-1480. 

1690. 

620.9 

A£/eq 

5.78 

0.77 

0.0000 

-1.51 

1.26 

H2O Adeq 

-0.08 

- 0.02 

0.002 

-0.04 

-0.03 

Ao; 

2207. 

-75.3 

-194.2 

-345.3 

398.0 

A_£/eq 

13.12 

4.25 

0.65 

0.001 

4.50 

HF Adeq 

-0.14 

-0.07 

- 0.02 

-0.0006 

-0.06 

Ao; 

3501. 

2272. 

1074. 

-247.5 

1650. 



FIG. 8 . d\P\{x) is calculated by finite difference d\P\{x) ~ 
, dxPx{x), of HF-^HzO at (a) d = 0.5 A for 
compression and at (b) d = 1.5 A for extension are plotted as 
a function of x and A. Nuclear positions are highlighted by red 
dashed lines, with F— >0 at a; = 0 A, H—^-void at a; = 0.93 A, 
two void —>-11 at X = —0.22 A and x = d, where void denotes 
the nuclei with zero charge. 


insets of Fig. [^b) and (c). The poor predictivey power 
of IPA has also recently been pointed out by Pulay and 
co-worker^^. By contrast, yields a very reason¬ 

able binding potential, albeit still far from being chemi- 

/o\ 

cally accurate. The superior performance of Aif^p, with 
respect to AEjpj^, indicates that the contributions of 
Coulomb and xc energy due to density response are cru¬ 
cial. In other words, matrix elements Jiajb and 'Kiajb in 
Eq. Q should not be neglected for non-vertical alchem¬ 
ical perturbations. 

Different predictive accuracy is found for compressing 
bonds d < deq versus stretching bonds d > deq- AE^p 
performs better in the region 0.5 A < d < 1.5 A. Simi¬ 
lar behaviour is also observed for other alchemical paths 
of compressing vs stretching bond. Also in this case, 
the aforementioned non-commutative asymmetric behav¬ 
ior of the predictions is observed. Namely, the AE^p 
based prediction for CH 4 —>■ NH 3 in Fig. IWc) is more 
accurate than for NH 3 —)■ CH 4 in Fig. [^bb Note that 
abrupt cha nges in electronic structure, as observed for 




Sec. 


IIIDl 


_^ ^ are not present when coupling these 

systems Since the accuracy of the second order es¬ 
timate is determined by how linearly the electron den¬ 
sity rearranges as a function of A, one expects a near¬ 
constant d\p for negligible higher order contributions. 
This is conhrmed through inspection of the integrated 
density response of the alchemical path HF—>H 2 0 in 
Fig.§ dxPxi x) varies less when A changes from zero to 
one for d = O.sA Fig.j^a), when compared with d = l.sA 
Fig. I^b). A near constant d\P\{x) at d = 0.5 A results 
in improved predictive accuracy. 


first order estimate does not even capture the changes in 
equilibrium bond length (Fig. [^b) and (cB. 

APjp^ yields a saddle point in Fig. Wa) and (d), 
instead of a minimum at optimized geometry. When 
the chemical compositions of Er and Ht are different, 

/q\ 

Aifip(^ results in in dramatic errors (worse than first or¬ 
der erstimates), as shown in the energy zoom out in the 


4- Integrated errors 

Table. HTT] summarizes the results for all 4x4 combi¬ 
nations of Hr —>■ idx, where mutual predictions of co¬ 
valent bond potentials in idx: {CH 4 , NH 3 , H 2 O, HF} 
are obtained based on only the single point wavefunc- 
tions obtained for the relaxed geometry of Hr-. {CH 4 , 
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NHa, H 2 O, HF}, respectively. This coupling matrix in 
chemical space (Table. Ill) is not symmetric (due to 
the non-commutative properties discussed above). Off- 
diagonal elements correspond to coupling paths involving 
changes in chemical composition and geometry. Diago¬ 
nal elements correspond to coupling paths that involve 
only changes in geometry, i.e. for the same stoichiome¬ 
try. Note that all error measures have been obtained via 
cubic spline fits. Therefore, also the predicted AE and 
the location of the energy minimum can be slightly non¬ 
zero even for the diagonal elements. These values should 
be considered noise: For the diagonal elements only the 
harmonic frequencies are meaningful. 

Table. HU confirms the trends observed above for first 
and second order. In general, best predictive power is 
found when the chemical composition of H^i is the same 
as i?T (diagonal elements). When the chemical com¬ 
position of Ht differs from the predictive accuracy 
deteriorates. This is not surprising and due to the per¬ 
turbing Coulomb potential being placed on the heavy 
atom in order to mutate it, e.g. from carbon to fluorine. 
Because of the strong accumulation of electron density 
(cusps) at the heavy atom’s site (6 to 7 electrons for car¬ 
bon to fluorine, respectively), this perturbation is quite 
severe. In the case of the diagonal element, by contrast, 
only the hydrogen atom is being annihilated and created, 
implying that the perturbing potential acts on the hydro¬ 
gen atom’s electronic density which is built up by only 
1 electron. This implies a less severe perturbation, and 
therefore worse predictive power can be expected for off- 
diagonal elements. 

The crucial importance of Coulomb and xc energy con¬ 
tribution to density response for second order alchemical 
perturbation is also confirmed for the other cases in Ta¬ 
ble. |III| These results clearly underscore the observation 
that IPA is a (very) poor approximation when it comes to 
estimate alchemical changes, yielding even worse predic¬ 
tions than the first order estimates. Interestingly enough, 
the first order estimate is even competitive in compari¬ 
son to the second order CP predictions. For example, 
using CH 4 as a reference compound the first order pre¬ 
diction deviates on average by -8.54 kcal/mol in the en¬ 
ergy, while AA^p deviates -15.99 kcal/mol. However, as 
the reference compound moves to the right hand side of 
the periodic table, the second order CP based estimate 
becomes more accurate than the first order based esti¬ 
mate. 


An additional aspect can be confirmed from inspec¬ 
tion of Table. |III[ The larger the perturbing potential, 
the worse the predictive accuracy of derivative based es¬ 
timate. More specifically, the larger the integrated norm 
of the difference between reference and target potential 
in the electronic Hamiltonian, the worse the predictive 
power. For example, using CH 4 as a reference compound, 
the prediction will be increasingly worse in the order of 
the respective predictions for NH 3 , H 2 O, and HF. Con¬ 
versely, using HF as a reference compound, the prediction 
will be increasingly worse in the order of the respective 


predictions for H 2 O, NH 3 , and CH 4 . Note that this is 
true for all first as well as second order estimates. 


IV. CONCLUSIONS 

The performance of truncated Taylor series for pre¬ 
dicting alchemical vertical changes in covalent bonding 
has been investigated in iso-electronic chemical spaces 
spanned by the external potentials of small molecules. 
For vertical linear transmutations (same geometry, same 
number of atoms, differing nuclear charges) our results 
suggest that chemical accuracy is possible when interpo¬ 
lating molecules containing p-block atoms from the third 
and fourth row using first order (Hellmann-Feynman the¬ 
ory) based predictions. Since first order estimates are an¬ 
alytical, this finding implies that one can scan potential 
energy surfaces of very many molecules with unprece¬ 
dented accuracy and speed as long as their stoichiome¬ 
tries are restricted to third and fourth row main group 
chemistries. First order based predictions of chemistries 
involving second row elements are only correct to a de¬ 
gree considered qualitative. 

Overall, we have found second order estimates to not 
provide sufficient improvement with respect to first or¬ 
der predictions (often even worse results) to warrant the 
investment in the additional overhead incurred. First 
order estimates are more accurate not because higher 
order terms are negligible, but rather due to the fact 
that (a) changes in relative energies (bonding) are al¬ 
ready near-linear (by optimizing the reference geometry) 
with respect to alchemical coupling (effectively canceling 
higher order terms), and (b) inflection points can occur 
which lead to worse predictions for second order esti¬ 
mates. For the interpolation of the pseudpotentials used 
in this study, inflection points near A = 0 are always ob¬ 
served when a lighter main group element is coupled to a 
heavier one. The absence of inflection points near A = 1 
improves the predictive power of the second order correc¬ 
tion: As such, the asymmetry of AA((i, dp) with respect 
to A = 0.5 results in asymmetric predictive performance. 

The choice of the reference geometry has a dramatic 
impact on the predictive power of the alchemical esti¬ 
mates. For covalent bond potentials, a linear relation¬ 
ship has been identified, (dopt ~ 0.76 djq -I- 0.97 A), that 
can be used to predict optimal dp requiring only rough 
estimates of the equilibrium bond-length in the target 
molecule (which can easily be obtained using universal 
force-fields or semi-empirical quantum chemistry meth¬ 
ods). 

We have found oscillating behaviour in the predictions 
of truncated Taylor series when varying the order in the 
non-vertical alchemical stretching of H^. The crossing of 
eigenvalue surfaces is due to the electron density’s neces¬ 
sity to be symmetric at A = 0 and A = 1. This leads to 
a diverging series Taylor series, yet the second order cor¬ 
rection could still provide fair predictions. The behav¬ 
ior of first and second order truncated alchemical Tay- 
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lor series expansions in non-vertical transmutations in 
chemical space has also been analyzed for molecules with 

ten electrons. Numerical evidence of the superior perfor- 

( 2 ) ( 2 ) 

mance of over AE^pj^ suggests that the response 

of Coulomb and xc energy to alchemical perturbation is 
crucial. 

In summary, our findings indicate that a careful 
choice of alchemical interpolation paths enables alchemi¬ 
cal derivatives to achieve predictive power with chemical 
accuracy for covalent bond potentials. Future work will 
deal with angles and torsions in lager molecules, as well 
as with solid metals and ionic crystals. 
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